First-principles calculations of phase transition, low elastic 
modulus, and superconductivity for zirconium 

Bao-Tian Wang/'^'Lj Peng Zhang,^ Han-Yu Liu,"^ Wei-Dong Li/ and Ping Zhang^'^'HI 

^Institute of Theoretical Physics and Department of Physics, 
—^ ■ Shanxi University, Taiyuan 030006, People's Republic of China 

^LCP, Institute of Applied Physics and Computational Mathematics, 
P.O. Box 8009, Beijing 100088, People's Republic of China 
~R • ^Department of Nuclear Science and Technology, 

Xi'an Jiaotong University, Xi'an 710049, People's Republic of China 
^ ' -^ State Key Lab of Superhard Materials, Jilin University, 

Changchun 130012, People's Republic of China 
^Center for Applied Physics and Technology, 
d • Peking University, Beijing 100871, People's Republic of China 

_L : Abstract 

Q ■ The elasticity, dynamic properties, and superconductivity of a, w, and /? Zr are investigated 

by using first-principles methods. Our calculated elastic constants, elastic moduli, and Debye 
temperatures of a and uj phases are in excellent agreement with experiments. Electron-phonon 
coupling constant A and electronic density of states at the Fermi level N{Ep) are found to increase 
Ti^lj- ■ with pressure for these two hexagonal structures. For cubic /3 phase, the critical pressure for 
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mechanical stability is predicted to be 3.13 GPa and at P=4 GPa the low elastic modulus (£'=31.97 
GPa) can be obtained. Besides, the critical pressure for dynamic stability of /3 phase is achieved 
by phonon dispersion calculations to be ~26 GPa. Over this pressure, A and N{Ep) of (3 phase 
decrease upon further compression. Our calculations show that the large value of superconducting 
transition temperature Tc at 30 GPa for /3 Zr is mainly due to the TAl soft mode. Under further 
compression, the soft vibrational mode will gradually fade away. 
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I. INTRODUCTION 



Group IV transition metals such as titanium, zirconium, and hafnium have wide apph 
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cations in the aerospace, nuclear, and biomedical industry [l|-|3|. The body-centered cubic 
(bcc) type (/3 phase) of Ti and Zr alloys in usage as metallic biomaterial have been studied 
extensively because of their properties of nontoxicity, high strength, good biocompatibility, 
and low elastic modulus [l|, \^. However, the investigations of the single crystal /3 Zr as 
biomaterial has not been reported. 

Actually, at ambient condition Zr crystallizes in hexagonal closed-packed (hep) structure 
{a phase). And at high temperature of 1135 K, it transforms into the /3 phase. At room 
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temperature and under pressure range of 2-7 GPa [5H7[ , the a phase transforms to another 
more open hexagonal structure of u phase with space group P6/mmm (No. 191). And under 
further high pressure of 30-35 GPa, the u^(3 phase transition has been observed [g, IsMll]. 
The Young's moduli of both a and u phases at ambient condition are largely bigger than 
that of a human bone (about 30 GPa). So the condition of low elastic modulus cannot be 
satisfied for these two phases. As for (3 phase, it cannot be exploited because the reverse 
phase transformations take place upon unloading. So it seems that we should not consider 
metal Zr as candidate biomaterial. But, fortunately recent experiment has successfully 
stabilized the /3 Zr at room temperature using compression stresses higher than 3 GPa 
plus the controlled application of shear |3|. This makes stabilizing the high temperature 
or pressure phases of Zr at ambient conditions reality. Therefore, the study of the elastic 
modulus for (3 Zr at low pressure is valuable. 

As another key objective of our work, we calculate the superconducting transition tem- 
perature Tc for metal Zr. Although the Tc of Zr has been observed in experiments |9(], 
to date no theoretical work is found in the literature. Our calculations will show that the 
increase (decrease) in T^ with pressure can be understood in terms of the enhancement (re- 
duction) of the electronic density of states (DOS) at the Fermi level and the corresponding 
behaviors of electron-phonon coupling constant. In present study, the dynamic stability of 
P Zr is tested carefully by performing phonon dispersion calculations. The /3— )-a; and (3^-a 
transitions are discussed in detail. The rest of this paper is arranged as follows. In Sec. II 
the computational methods are described. In Sec. Ill we present and discuss our results. In 
Sec. IV we summarize the conclusions of this work. 



II. COMPUTATIONAL METHODS 



A. Computational details 



The first-principles density functional theory (DFT) calculations on the basis of the 



frozen-core projected augmented wave (PAW) method of Blochl 12| are performed within 



the Vienna ab initio simulation package (VASP) |13i], where the Perdew, Burke, and Ernz- 
erhof (PEE) [ij] form of the generalized gradient approximation (GGA) is employed to 
describe electron exchange and correlation. For the plane- wave set, a cutoff energy of 500 
eV is used. The P-centered k point-meshes in the full wedge of the Brillouin zone (BZ) are 
sampled by 18x18x16, 16x16x18, and 18x18x18 grids according to the Monkhorst-Pack 
(MP) [l5| for a (two-atom cell), u (three-atom cell), and (3 (two-atom cell) Zr, respectively, 
and all atoms are fully relaxed until the Hellmann-Feynman forces become less than 0.001 
eV/A. The Zr Ad^bs^ orbitals are explicitly included as valence electrons. The pseudopoten- 
;ial plane- wave method within the PBE-GGA through the QUANTUM-ESPRESSO package 
16| is employed to calculate the electronic properties, lattice dynamics, and electron-phonon 
coupling (EPC) for Zr, where the Zr 4s^4p^4(i^5s^ are treated as valence electrons. Conver- 
gence tests give the choice of kinetic energy cutoffs of 60 Ry with Gaussians width of 0.05 
Ry for all three phases of Zr and 12x12x8, 8x8x12, and 16x16x16 MP grids of A:-point 
meshes for a, u, and /3 Zr, respectively. Phonon frequencies are calculated based on the 
density functional linear- response method [17|,ll8|. In the interpolation of the force constants 
for the phonon dispersion curve calculations, 6x6x4, 4x4x6, and 8x8x8 g-point meshes 
and denser 24x24x16, 16x16x24, and 32x32x32 A:-point meshes in the first BZ are used 
for a, oj, and (3 phases of Zr, respectively. 

B. Mechanical properties 

To avoid the Pulay stress problem, the geometry optimization at each volume is per- 
formed at fixed volume rather than constant pressure. The theoretical equilibrium volume, 
bulk modulus B, and pressure derivative of the bulk modulus B' are obtained by fitting 
the energy- volume data in the third-order Birch- Murnaghan equation of state (EOS) |l9 |. 
Elastic constants for cubic symmetry (Cn, Cu, and C44) and hexagonal structure (Cn, 
C12, Ci3, C33, and C44) are calculated by applying stress tensors with various small strains 



onto the equilibrium structures. The strain amphtude S is varied in steps of 0.006 from 
(5=— 0.036 to 0.036. After obtaining elastic constants, the polycrystalline bulk modulus B 
and shear modulus G are calculated from the Voigt-Reuss-Hill (VRH) approximations J20 |. 
The Young's modulus E and Poisson's ratio v are calculated through E = 9BG/{3B + G) 
and V = {3B — 2G)/[2{3B + G)]. In calculation of the Debye temperature {9d), we use the 
relation 
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where h and ks are Planck and Boltzmann constants, respectively, n is the number of atoms 
in the molecule, Q is molecular volume, and Vm is the average sound wave velocity. The 
average wave velocity in the polycrystalline materials is approximately given as 
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where Vt = \/G j p (p is the density) and vi = ^/{3B + AG)/3p are the transverse and 
longitudinal elastic wave velocity of the polycrystalline materials, respectively. 



C. Superconduntivity 



The superconducting transition temperature Tc is evaluated by using the Allen-Dynes 
modified McMillan equation 2M, l22 | 
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is the logarithmic average frequency. 
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is the electron-phonon coupling constant, and p* is the Coulomb pseudopotential. The 
Eliashberg electron-phonon spectral function a^E{(jj) is written as 
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where N{Ep) is the electronic density of state (DOS) at the Fermi level. 
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FIG. 1: (Color online) Calculated enthalpy differences of a and w phases with respect to /3 phase 
as functions of pressure. 

III. RESULTS 



At K, the Gibbs free energy is equal to the enthalpy H. After calculation, we can plot 
in Fig. [1] the relative enthalpies of the a and the u phases with respect to the (3 phase as 
functions of pressure. The crossing between the a and u enthalpy curves readily gives phase 
transition pressure of —3.7 GPa, which indicates that at ambient pressure the u phase is 
more stable than the a phase. This fact is in d isag reement with experiment, but coincides 



well with recent DFT-PBE [23| and FP-LAPW 2J| calculations. The disagreement between 
theory and experiment mainly originates from temperature effect. Actually, previous DFT- 
PBE 23J study clearly shows that the transition of a^u occurs at 1.7 GPa at r=300 K. 



As for the uj^(3 transition, a pressure of 32.4 GPa is obtained. This value falls exactly in 
the pressure range of 30—35 GPa measured by experiments |6|, IsMllj and also fairly accords 



with other calculations 
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In Table HI we report our calculated lattice parameters, bulk modulus, pressure derivative 



of the bulk modulus, and elastic constants o 



comparison, the experimental values 
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a, u, and B Zr at different pressures. 

n fin 

27l | and other theoretical results [l|, |28|, 
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are also listed. Obviously, our calculated results accord well with experiments and other 
calculations for a and u phases. However, for /3 Zr the difference of lattice parameter and 
bulk modulus between our calculation and the experiment from Zhao et al. 7] is evident. 
This is due to the fact that the experiment is conducted at r=973 K, while our calculation 



TABLE I: Calculated lattice parameters (a or c), bulk modulus (B), pressure derivative of the 
bulk modulus {B ), and elastic constants of a, uj, and /3 Zr at different pressures. For comparison, 
experimental values and other theoretical results are also listed. 
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This work 





3.236 


5.168 


96.0 


3.12 
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71.0 
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26.0 






5.35 
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FIG. 2: (Color online) Calculated compression curves of /3 Zr compared with experimental mea- 
surements. 



is valid only at K. Similar temperature effect for P Zr is also found in studying the 
compression behaviour. As shown in Fig. [21 although good consistence with the experiment 
[id ] at high pressure domain is clear, the values of the atomic volumes are smaller than 
that by experiment at T=973 K. Nevertheless, the compression performance of /3 Zr is 
wholly consistent with experiments. From Table [U one can find that the difference of elastic 
constants at zero pressure between calculation and experiment is evident. This also can be 
attributed to the temperature effect. In fact, the good accordance of the lattice parameter 
and elastic constants of /3 Zr at GPa between our calculation and previous DFT-PBE [l| 
study supplies the safeguard for our study of structure, mechanical, and electronic properties 
of P Zr. Additionally, the mechanical stability of a and cj Zr at GPa and at some typical 
finite pressures can be predicted from the elastic constants data. But the elastic constants 
of /3 Zr at GPa illustrate that the cubic phase is mechanically unstable, which is in good 
agreement with the results by Ahuja et al. 29| and Ikehata et al. [Ij. Actually, the 
mechanical stability of /3 Zr can be obtained through applying external pressure and our 
data show that along with the increase of pressure from GPa to 60 GPa, the value of 
Cii — Ci2 increases almost linearly from —5.7 GPa to 157 GPa. Fitting the curve (not 
shown) of the pressure behaviour of Cu — C12 by first-order polynomial function, we find 
that the value of Cu — C12 becomes positive at P >3.13 GPa. In fact, at P=4 GPa our 
first-principles calculated Cu, C12, and C44 equal to 101.8, 98.2, and 28.7 GPa, respectively, 
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TABLE II: Calculated elastic moduli, Poisson's ratio {v), density (p), transverse {vt), longitudinal 

(f/), and average {vm) sound velocities calculated from poly crystalline bulk and shear modulus, 
and Debye temperature of a, w, and /3 Zr at different pressures. For comparison, experimental 
values and other theoretical results are also listed. 

Phase Method Pressure B G E v p vt vi Vm Gd 

(GPa) (GPa) (GPa) (GPa) (g/cm^) (km/s) (km/s) (km/s) (K) 

97.4 33.8 91.0 0.344 6.465 2.288 4.695 2.571 267.3 

110.4 34.9 94.8 0.357 6.804 2.266 4.803 2.550 269.8 
95.3 36.3 96.6 0.331 6.511 2.361 4.698 275.9 

97.5 43.6 113.7 0.306 6.568 2.575 4.868 2.878 300.9 

113.5 49.0 128.5 0.311 6.948 2.655 5.074 2.970 276.4 
10.01 123.2 51.6 135.8 0.316 7.177 2.681 5.172 3.001 282.3 

Expt.'' 104.0 45.1 118.3 0.311 6.589 2.616 4.996 306.1 

6.88 2.638 5.177 

10.59 2.652 5.284 

/3 This work 10 112.1 18.7 53.1 0.421 7.309 1.599 4.330 1.816 196.7 

20 131.2 29.8 83.2 0.394 7.889 1.944 4.655 2.199 244.4 

30 148.4 39.8 109.6 0.377 8.415 2.175 4.893 2.454 278.7 

35 156.5 45.5 124.6 0.367 8.662 2.293 5.008 2.584 296.3 

40 164.0 50.9 138.3 0.360 8.900 2.390 5.104 2.691 311.4 

50 179.3 58.3 157.8 0.353 9.353 2.497 5.243 2.809 330.4 

60 193.6 68.0 182.5 0.343 9.781 2.636 5.391 2.961 353.5 

Expt.* 96.7 18.1 51.2 0.412 6.420 1.681 4.339 1.906 197.8, 177'= 
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"' Reference [30], ^ Calculated using present scheme with elastic constants at T=915°C from 
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reference 12611. '= Reference 
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which exphcitly indicate the elastic stabihty of bcc Zr under this pressure. 

After obtaining elastic constants, the elastic moduli, Poisson's ratio (v), transverse sound 
velocities (vt), longitudinal sound velocities (vi), average sound velocities (vm), and Debye 
temperature of a, u, and Zr at different pressures are calculated and tabulated in Table [ITl 



The experimental values 26|, l30|] are also collected in Table HTl for comparison. The excellent 
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FIG. 3: (Color online) Phonon dispersion at GPa (a) and phonon dispersion, Eliashberg phonon 
spectral function, and phonon DOS at 26, 30, and 60 GPa (b) for /3 Zr. 

coincidence between our calculation and experiment for a and u phases can be seen. For /3 
phase, the bulk modulus 5, shear modulus G, Young's modulus E, and Debye temperature 
all linearly increase with pressure. The Young's modulus E at P=4 GPa (just over the 
critical pressure 3.13 GPa for the elastic stability of /3 Zr) is calculated to equal to 31.97 
GPa, which is very near the Young's modulus value (~30 GPa) of a human bone. Thus, we 
successfully indicate the importance of stabilizing the high temperature or pressure phase 
of Zr at ambient conditions. In Table [TTl we also list the results of /3 phase calculated using 
present scheme with elastic and lattice constants at T=915°C from Ref. 26|]. The derived 
Debye temperature is slightly higher than that given in Ref. [26|. 



Phonon dispersion results 
the /3 phase of Ti and Hf 



31 



GiB Zr at selected pressures are shown in Fig. 121 Like 



32| . the high pressure /3 phase of Zr is also dynamically 



unstable at low pressure [Fig. [3]^a)]. Our phonon dispersion calculations at high pressure 
illustrate that the /3 Zr becomes dynamically stable at around 26 GPa. In fact, no imaginary 
frequencies are observed at 26-60 GPa pressure domain as clearly indicated in Fig. [3]^b). 
Therefore, although recent experiment by Perez-Prado and Zhilyaev [3] has successfully 
stabilized the high pressure/temperature phase of Zr at ambient condition, our phonon 
dispersion calculations illustrate that the /3 phase is dynamically unstable in pressure domain 
of 0-25 GPa at K. This fact confirms that the temperature effect on the phase stability is 
also critical for the experimentally observed J3| /3 phase at ambient condition. Besides, the 
dynamical stability test at 26-60 GPa supplies the safeguard for our following study of the 
superconductivity of /3 Zr at high pressure. As shown in Fig. [3]^a), our phonon dispersion 
results indicate soft modes for the longitudinal L |[111] and transverse T |[110] phonons. 
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FIG. 4: (Color online) Schematic illustrations of the structural transition for (a) f3^>-uj and (b) 



As indicated by experiment [26|, the L 3 [111] unstable phonon branch is responsible for 
the P^uj transition and the T |[110] implies the P—^a transition. In Fig. |U we show the 
paths of P^uj and P^^a transitions. As shown in Fig. Hl^a), the bcc /3 phase can be viewed 
as ABC periodically layered structure along the [HI] direction and there are two periods 
of ABC layers in its unit cell. For u phase, it can be viewed as AB periodical layered 
structure along the [0001] direction [Fig. Hl^a)]. The atoms on B and C layers in /3 phase 
move easily along the soft mode direction of [HI] and can transfer to be the B layer of cu 
phase. Therefore, the atomic displacement along the [111] direction is the main origin for 
the /3— >a; phase transition. For /3— >a; transition, one can regard the (3 phase as two layered 
(AB) structure along the [llO] orientation [Fig. lU^b)]. In transition, firstly the atoms on 
each layer expand along the [110] direction to form a hexagon, and second, the adjacent 
(110) planes slip relatively along the [110] direction to create the hep structure. 

For high pressure phase of Zr, it has evident predominance for high-temperature super- 
conductivity. The EPC calculations have been performed to explore the superconductivity. 
The Eliashberg phonon spectral function a'^F{uj) and phonon DOS at 26, 30, and 60 GPa 
for /3 Zr are shown in Fig. |3]^b). Logarithmic average phonon frequency wiog and the EPC 
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TABLE III: Calculated logarithmic average of vibrational frequencies a;iog and EPC constant A of 
a, CO, and /3 Zr at different pressures. 

Phase Pressure ujiog A 

(GPa^ m 

a 

5 
uj 

10 

20 

30 
P 26 

28 

30 

40 

50 

60 

strength A of a, a; , and (3 Zr at different pressures are listed in Table ITTTl The superconduct- 
ing Tc has been estimated from Eqs. (3), and a typical value of Coulomb pseudopotential 
fi*=0.12 is used. Our calculated Tc, together with N{Ep), as a function of pressure for the 
three phases are plotted in Fig. [5l Note that there are two different kinds of atoms for 
u Zr: Zrl in la and Zr2 in 2d sites. Clearly, the predicted T^ is in good agreement with 
experimental results. As indicated by McMilan equation [Eqs. (3)], Tc has tight relationship 
with Uiog, A, and N{Ep). From 0-7 GPa, total N{Ep) of a phase increases with pressure. In 
pressure range of 7-30 GPa, the average N{Ep) of u phase over two kinds of atomic arrange- 
ments also exhibits increasing behavior upon compression. In contrast, total N{Ep) of /3 Zr 
decreases with pressure in the whole domain of 26-60 GPa. All these pressure-dependent 
behaviors of N{Ep) obey one known fact that a more stable phase has a lower value of 
N{E-p). Through plotting partial N{Ep) of s, p, and d orbitals (not shown), we find that 
the main contribution to the pressure-dependent behavior of the total N{Ep) comes from 
d orbital. The s orbital has no contribution and the contribution from p orbital is limited. 
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FIG. 5: (Color online) Calculated total density of state at the Fermi level N{Ep) and supercon- 
ducting transition temperature Tc as functions of pressure for a, u, and /? Zr. The experimental 
values of Tc from Ref. [3] are also presented. The two shaded resigns indicate the a— )-a; and w— )-/S 
phase transition ranges of 2—7 GPa J5|47|] and 30—35 GPa 
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From Fig. [5l one can see that the experimental Tc changes in almost the same way as that 
of the N{Ep). So, it seems that relative larger value of the N{Ep) for /3 Zr at 26 GPa may 
predict larger T^ for superconductivity. However, N{Ep) is less directly important than uj\og 
and A for T^. 

From Table llllt one can see that the calculated A for a and u phases increases with pres- 
sure while uiog does not show any clear dependence on pressure. As a result, the increasing 
behavior of Tc upon compression is mainly due to the increase of A. For /3 Zr, wiog and A 
have an almost opposite evolution with pressure. However, the changing rate for u\og and 
A is different. The increasing rate of u\og is clearly larger than the decreasing rate of A in 
pressure range of 26 to 30 GPa. After being compressed over 30 GPa, the decreasing rate 
of A exceeds the increasing rate of wiog- So the main affection to the pressure performance 
of Tc comes from the increase of wiog in pressure domain 26-30 GPa, while in 30-60 GPa 
is from the decrease of A. Therefore, the relatively small value of Tc at 26 GPa for /3 Zr 
is understandable. Furthermore, the lowest transverse-acoustics (TAl) soft mode for bcc /3 
Zr in pressure range of 26-30 GPa is evident [Fig. [3](b)] and can lead a significant electron- 
phonon contribution. Therefore, as indicated in Table UTTl the EPC constant A at 26-30 GPa 
is predominantly larger than that of both a and u phases. And, along with increasing the 
pressure, soft vibrational mode for (3 Zr gradually fade away, resulting in the decrease of A 



12 



and Tc. Thus, the decrease of A and Tc for /3 Zr under pressure mainly originates from the 
fading of soft vibrational modes induced by increasing pressure. 

IV. CONCLUSION 

In summary, the elasticity, dynamic properties, and superconductivity of Zr under pres- 
sure up to 60 GPa have been studied by means of the first-principles DFT-PBE method. 
Our results have shown that the structural parameters, elastic constants, elastic moduli, 
Poisson's ratio, ultrasonic velocities, and Debye temperature of a and u phases coincide 
well with experiments. The a^u and w— t-/? phase transition pressures at T=0 K have been 
calculated to be —3.7 GPa and 32.4 GPa, respectively. The compression curve of f3 Zr is 
well consistent with experiments. Our calculations have explicitly indicated that the low 
elastic modulus of Zr occurs at low pressure in /3 phase. 

In phonon dispersion study, the soft modes of f3 phase along [111] and [110] directions 
have been shown and the /3— ^w and /3— >q; transition paths have been predicted. Dynamic 
stability test has illustrated that f3 Zr is unstable till being compressed over 25 GPa. In 
addition, superconductivity has been obtained by electron-phonon coupling calculations and 
our calculated T^ accords well with experiments. Through analyzing the electronic density 
of states at the Fermi level, we have derived that the main contribution to the pressure- 
dependent behavior of the superconductivity comes from the d orbital. The large T^ at 
around 30 GPa for /3 Zr is mainly due to the TAl soft mode. Under pressure, the increase or 
decrease of Tc for Zr in all three phases has tight relation with the corresponding behavior 
of the EPC constant A. 
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